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ABSTRACT 

We explore the utility of Karhunen Loeve (KL) analysis in solving practical problems in the analysis 
of gravitational shear surveys. Shear catalogs from large-field weak lensing surveys will be subject to 
many systematic limitations, notably incomplete coverage and pixel-level masking due to foreground 
sources. We develop a method to use two dimensional KL eigenmodes of shear to interpolate noisy 
shear measurements across masked regions. We explore the results of this method with simulated 
shear catalogs, using statistics of high-convergence regions in the resulting map. We find that the 
KL procedure not only minimizes the bias due to masked regions in the field, it also reduces spurious 
peak counts from shape noise by a factor of ~ 3 in the cosmologically sensitive regime. This indicates 
that KL reconstructions of masked shear are not only useful for creating robust convergence maps 
from masked shear catalogs, but also offer promise of improved parameter constraints within studies 
of shear peak statistics. 

Subject headings: gravitational lensing — dark matter — large-scale structure of universe 



1. INTRODUCTION 

Since its discovery over a decade ago, cosmic 
shear - the coherent gravitational distortion of light 
from distant galaxies due to large scale structure - 
has beco me an important tool in pre cision cosmol- 
ogy, (see iBartelmann fe Schneider! 12001. for a review). 
Kaiser et al.l (|1995l) first proposed a method of using 
weak gravitational lensing to probe the projected gravi- 
tational mass within an observed field. This method and 
its variations have often been applied to reconstructions 
of clusters using deep observations within small (a few 
square arcminutc) fields. In this situation, the measured 
shear often dominates the shape noise, and useful priors 
can be applied based on assumptions regarding the lo- 
cation of the cluster center, symmetries of cluster mass 
profile, and the correlation of mass with light. 

Currently, a new generation of wide-field weak lens- 
ing surveys are in the planning and construction stages. 
Among them are the the Dark Energy Survey (DES), 
the Panoramic Survey Telescope & Rapid Response Sys- 
tem (PanSTARRS), the Wide Field Infrared Survey Tele- 
scope (WFIRST), and the Large Synoptic Survey Tele- 
scope (LSST), to name a few. These surveys, though 
not as deep as small- field space-based lensing surveys, 
will cover orders-of-magnitude more area on the sky: up 
to ~ 20, 000 square degrees in the case of LSST. 

This is a fundamentally different regime than early 
weak lensing reconstructions of single massive clusters: 
the strength of the shear signal is only ~ 1%, and is dom- 
inated by ~ 30% intrinsic shape noise. This, combined 
with source galaxy densities of only n ~ 20 — 50 arcmin -2 
(compared with n > 100 arcmin -2 for deep, space- 
based surveys) and atmospheric psf effects leads to a 
situation where the signal is very small compared to 
the noise. Additionally, in the wide-field regime, the 
above-mentioned priors cannot be used. Nevertheless, 



many methods have been developed to extract useful 
information from wide-field cosmic shear surveys, in- 
cluding measurin g the N-point power spectra and corre- 
lation functions dSchneider et al.l [2002at iTakada fc JainI 
120041 : iHikage et al.l 120101 ), performing log tran sforms 
of the convergence field (iNevrinck et alJ 120091 120101 : 
iScherrer et aljf20l"ot iSeo et al.ll201ll ). analy zing statistics 
of convergence and aperture m a ss peaks dMarian et al.l 
2Ol0l [Dietoich fc Hartlapl [2011 ISchmidt fc Rozol 120101: 



Kratochvil et al.l I201fl iMaturi et al.l 120111 ). Another 
well-motivated application of wide-field weak lensing is 
using wide-field mass reconstructions to minimize the ef- 
fect mass-sheet degeneracy in halo mass determination. 

Many of the above applications require reliable recov- 
ery of the projected density, either in the form of the 
convergenc e k, or filter-based q uantities such as aper- 
ture mass (jSchneider et al.lll998l ). Because each of these 
amounts to a non-local filtering of the shear, the pres- 
ence of masked regions can lead to a bias across sig- 
nificant portions of the resulting maps. Many of these 
methods have been demonstrated only within the context 
of idealized surveys, with exploration of the complica- 
tions of real-world survey geometry left for future study. 
Correction for masked pixels has been studied within 
the c ontext of shear po wer spectra (jSchneider et alJ 
120101 : IHikage et al.l 120101) but has not yet been sys- 
tematically addressed within the context of mapmak- 
ing and the associated stat i stical methods (se e, however, 
iPadmanabhan et all I2001 : iPires et alJ [200l for some 
possible approaches). We propose to address this missing 
data problem through Karhunen-Loeve (KL) analysis. 

In Section [2] we summarize the theory of KL analysis 
in the context of shear measurements, including the use 
of KL for interpolation across masked regions of the ob- 
served field. In Section [3] we show the shear eigenmodes 
for a particular choice of survey geometry, and use these 
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eigenmodcs to interpolate across an artificially masked 
region in a simulated shear catalog. In Section [4] we dis- 
cuss the nascent field of "shear peak statistics" , the study 
of the properties of projected density peaks, and propose 
this as a test of the possible bias imposed by KL analysis 
of shear. In Section [5] we utilize simulated shear cata- 
logs in order to test the effect of KL interpolation on the 
statistics of shear peaks. 

2. KARHUNEN-LOEVE ANALYSIS OF SHEAR 

KL analysis is a commonly used statistical tool in a 
broad range of astronomical applic ations, from, e.g. stud- 
ies of galaxy and quasa r spectra < Connolly et al.l 119951 : 
iConnollv fc Szalavlll99llYip et al Il2004allbl). to~analysis 
of th e spatial distribut ion of galaxies (iVogeley fc Szalavl 
[1991 IMatsubara et al.ll2000t iPope et al.l 12004) . to char- 
acterizat ion of the expected err o rs in weak lensing 
surve ys (jKilbinger fc Munsl ii 2006; Mun shi fc Kilbingerl 
120061 ). Any set of iV-dimensional data can be represented 
as a sum of TV orthogonal basis functions: this amounts 
to a rotation and scaling of the TV-dimensional coordi- 
nate axis spanning the space in which the data live. KL 
analysis seeks a set of orthonormal basis functions which 
can optimally represent the dataset. The sense in which 
the KL basis is optimal will be discussed below. For the 
current work, the data we wish to represent are the ob- 
served gravitational shear measurements across the sky. 
We will divide the survey area into N discrete cells, at lo- 
cations Xi, 1 < i < N. From the ellipticity of the galax- 
ies within each cell, we infer the observed shear j°(xi), 
which we assume to be a linear combination of the true 
underlying shear, j(xi) and the shape noise n 7 (ai;)[!] In 
general, the cells may be of any shape (even overlapping) 
and may also take into account the redshift of sources. 
In this analysis, the cells will be square pixels across the 
locally flat shear field, with no use of source redshift in- 
formation. For notational clarity, we will represent quan- 
tities with a vector notation, denoted by bold face: i.e. 



7 = [7i,72 



IT. 



7, = 7(034). 
2.1. KL Formalism 



KL analysis provides a framework such that our mea- 
surements 7 can be expanded in a set of N orthonormal 
basis functions {&j(xi), j = 1,-ZV}, via a vector of coef- 
ficients a. In matrix form, the relation can be written 



7 = *a 



(1) 



where the columns of the matrix \I r are the basis vectors 
if?i. Orthonormality is given by the condition = 
Sij , so that the coefficients can be determined by 



7 



(2) 

A KL decomposition is optimal in the sense that it seeks 
basis functions for which the coefficients arc statistically 
orthogonal^ that is, they satisfy 



( a i)k 



(3) 



1 Throughout this work, we assume we are in the regime where 
the convergence k < 1 so that the avera ge observed ellipticity in a 
cell is an unbiased estimator of shear; sec Bartclmann & Schneider 

2 Note that statistical orthogonality of coefficients is conceptu- 
ally distinct from th e geometric orthogonalit y of the basis func- 
tions themselves; sec Vogclcv & Szalay (1996) for a discussion of 
this property. 



where angled braces (• • •) denote averaging over all real- 
izations. This definition leads to s everal important prop- 
erties fsee lVogelev fc Szalay 1996, for a thorough discus- 
sion & derivation): 

1. KL as an Eigenvalue Problem: Defining the 
correlation matrix £y = (7^7*), it can be shown 
that the KL vectors are eigenvectors of £ with 
eigenvalues A; = (a 2 ). For clarity, we'll order the 
cigenbasis such that Xi > A; + i V i 6 (1, TV — 1). We 
define the diagonal matrix of eigenvalues A, such 
that Aij = XiSij and write the eigenvalue decom- 



position in compact form: 
£ = *A* f 



(4) 



2. KL as a Ranking of Signal-to-Noise It can 

be shown that KL vectors of a whitened covariance 
matrix (see Section |2~2|) diagonalize both the signal 
and the noise of the problem, with the signal-to- 
noise ratio proportional to the eigenvalue. This 
is why KL modes are often called "Signal-to-noise 
eigenmodcs" . 

3. KL as an Optimal Low-dimensional Rep- 
resentation: An important consequence of the 
signal-to-noise properties of KL modes is that the 
optimal rank-n representation of the data is con- 
tained in the KL vectors corresponding to the n 
largest eigenvalues: that is, 



n<N 



7 



00 



(5) 
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minimizes the reconstruction error between 7^ 
and 7 for reconstructions using n orthogonal ba- 
sis vectors. This is the theoretical basis of Prin- 
cipal Component Analysis (sometimes called Dis- 
crete KL), and leads to a common application of 
KL decomposition: filtration of noisy signals. For 
notational compactness, we will define the trun- 
cated eigenbasis and truncated vector of co- 
efficients a(„) such that Equation [5] can be written 

in matrix form: 7^™^ = iff^a^y 

2.2. KL in the Presence of Noise 

When noise is present in the data, the above proper- 
ties do not necessarily hold. To satisfy the statistical 
orthogonality of the KL coefficients a and the result- 
ing signal-to-noise properties of the KL eigenmodes, it 
is essential that the noise in the covariance matrix be 
"white": that is, A/" 7 = (n T nD cx I. This can be ac- 
complished through a judicious choice of binning, or by 
rescaling the covariance with a whitening transformation. 
We take the latter approach here. 

Defining the noise covariance matrix A/" 7 as above, 
the whitened covariance matrix can be written £ w = 
Af~ 1/2 £J\f~ 1/2 . Then the whitened KL modes become 

'I'vkAw^vi' = £w The coefficients a\y are calculated 
from the noise-weighted signal, that is 



a w = *| 4/ A/' 7 1/2 (7 + n 7 ) 



(0) 
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For the whitened KL modes, if signal and noise are un- 
correlated, this leads to (awo, w ) = Aw + I- that is, 
the coefficients aw are statistically orthogonal. For the 
remainder of this work, we will drop the subscript u w" 
and assume all quantities to be those associated with the 
whitened covariance. 

2.3. Computing the Shear Correlation Matrix 

The KL reconstruction of shear requires knowledge of 
the form of the pixel-to-pixel correlation matrix In 
ma ny applications of KL (e.g. analysis of galaxy spec- 
tra, IConnollv et"aT1[T995h this correlation matrix is de- 
termined empirically from many realizations of the data 
(i.e. the set of observed spectra). In the case of weak 
lensing shear, we generally don't have many realizations 
of the data, so this approach is not tenable. Instead, we 
compute this correlation matrix analytically. The corre- 
lation of the cosmic shear signal between two regions of 
the sky Ai and Aj is given by 



£ij = <7i7j) + ( n i n j) 

= / d 2 Xi I d 2 Xj^+(\xi — Xj 

J A, J A, 



Sii-t (7) 



where er e is the intrinsic shape noise (typically as- 
sumed to be ~ 0.3), n is the average galaxy count per 
pixel, and (9) is the "+" shear correlation function 
(ISchneider etall l2002ah . £+(0) can be expressed as an 
integral over the shear power spectrum: 



£f(0) = ^ J°° d£ £P^(£)J, 



{16) 



(8) 



where Jo is the zeroth-order Bessel function of the first 
kind. The shear power spectrum P 1 {€) can be expressed 
as an appropriately weighted line-of-sig ht integral over 
the 3 D mass power spectrum (see, e.g. iTakada fc Jainl 
2004): 



P 7 {£) 



dxW 2 (x)x- 2 P8[k = -;z( X ) 



(9) 



Here x is the comoving distance, Xs is the distance to 
the source, and W{x) 1S the lensing weight function, 

W(x) = ; v — / dzn{z) — — (10) 

2a(x) rig J z[x) x[z) 

where n(z) is the redshift distribution of galaxies. We as- 
sume a DES-like survey, where n(z) has the approximate 
form 

n(z) oc z 2 exp[-(z/z ) 1 - 5 ] (11) 

with zq = 0.5, where n(z) is normalized to the observed 
galaxy density n g = 20 arcmin - . 

The 3D mass power spectrum Pg(k,z) in Equation [5] 
can be computed theoretically. In this work we com pute 
Ps(k, z) using the halo model of lSmith et al.l (120031 ). and 
compute the correlation matrix £ using Equations Ffl llll 
When computing the double integral of Equation we 
calculate the integral in two separate regimes: for large 
separations (9 > 20 arcmin), we assume £+(9) doesn't 
change appreciably over the area of the pixels, so that 
only a single evaluation of the x+ {&) is necessary for each 



pixel pair. For smaller separations, this approximation 
is insufficient, and we evaluate ^ using a Monte-Carlo 
integration scheme. Having calculated the theoretical 
correlation matrix £ for a given field, we compute the 
KL basis directly using an eigenvalue decomposition. 

2.4. Which Shear Correlation? 

Above we note that the correlation matrix of the mea- 
sured shear can be expressed in terms of the "+" correla- 
tion function, £ + (9) . This is not the only option for mea- 
surem ent of shear correlations (see, e.g. ISchneider et al.l 
l2002af ). So why use £+(6) rather than £-(9)7 The answer 
lies in the KL formalism itself. The KL basis of a quan- 
tity 7 is constructed via its correlation (77^). Because of 
the complex conjugation involved in this expression, the 
only relevant correlation function for KL is £ + (9) by defi- 
nition. Nevertheless, one could object that by neglecting 
KL under- utilizes the theoretical information avail- 
able about the correlations of cosmic shear. However, in 
the absence of noise, the two correlation functions con- 
tain identical information: cither function can be deter- 
mined from the other. In this sense, the above KL for- 
malism uses all the shear correlation information that is 
available. 

One curious aspect of this formalism is that the theo- 
retical covariance matrix and associated eigenmodes are 
real- valued, while the shear we are trying to reconstruct 
is complex- valued. This can be traced to the computa- 
tion of the shear correlation: 

£+ = (77*) = (7t7t) + (7x7x)+i[(7t7x)-(7x7t)] (12) 

By symmetry, the imaginary part of this expression is 
zero. At first glance, this might seem a bit strange: how 
can a complex- valued data vector be reconstructed from 
a real-valued orthogonal basis? The answer lies in the 
complex KL coefficients a^: though each KL mode con- 
tributes only a single phase across the field (given by 
the phase of the associated ai), the reconstruction has 
a plurality of phases due to the varying magnitudes of 
the contributions at each pixel (given by the elements of 
each basis vector ^,). 

An important consequence of this observation is that 
the KL modes themselves are not sensitive by construc- 
tion to the E-modc (curl-free) and B-mode (divergence- 
free) components of the shear field. As we will show be- 
low, however, the signal-to-noise properties of KL modes 
lead to some degree of sensitivity to the E and B-mode 
information in a given shear field (See Section [5} . 

2.5. Interpolation using KL Modes 

Shear catalogs, in general, are an incomplete and in- 
homogeneous tracer of the underlying shear field, and 
some regions of the field may contain no shear informa- 
tion. This sparsity of data poses a problem, because the 
KL m odes are no longer o rthog onal over the incomplete 
field. IConnollv fe Szalavl (|1999f ) demonstrated how this 
missing-information problem can be addressed for KL de- 
compositions of galaxy spectra. We will summarize their 
results here. First we define the weight function w(xi). 
The weight function can be defined in one of two ways: a 
binary weighting convention where w(xi) = in masked 
pixels and 1 elsewhere, or a continuous weighting conven- 
tion where w(xi) scales inversely with the noise [A/" 7 ]m. 
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The binary weighting convention treats the noise is part 
of the data, and so the measurements should be whitened 
as outlined in Section [2T2J The continuous weighting con- 
vention assumes the noise is part of the mask, so data 
and noise are not whitened. We find that the two ap- 
proaches lead to qualitatively similar results, and choose 
to use the binary weighting convention for the simplicity 
of comparing masked and unmasked cases. 

Let 7 be the observed data vector, which is uncon- 
strained where w(xi) = 0. Then we can obtain the KL 
coefficients ai by minimizing the reconstruction error of 
the whitened data 

X 2 = (Af- 1/2 1° - ^^a^WiM- 1 ' 2 ^ - * (n) a ( „)) 

where we have defined the diagonal weight matrix Wij = 
w(xi)Sij. Minimizing Equation 1131 with respect to a 
leads to the optimal estimator a, which can be expressed 

aw^M^l^WAr-^Y (14) 

Where we have defined the mask convolution matrix 
7Vf(„j = V&t nW\&(„). These coefficients d(„) can then 
be used to construct an estimator for the unmasked shear 
field: 

7 W = A^y 2 * (n) a ( „ ) (15) 

In cases where the mask convolution matrix M< n ) is sin- 
gular or nearly singular, the estimator in Equation [T5] 
can contain unrealistically large values within the recon- 
struction j( n \ This can be addressed either by reducing 
n, or by adding a penalty function to the right side of 
Equation ll3l One convenient form of this pe nalty is the 
generalized Wiener filter (see lTegmarklll997l) . which pe- 
nalizes results which deviate from the expected correla- 
tion matrix. Because the correlation matrix has already 
been computed when determining the KL modes, this fil- 
ter requires very little extra computation. With Wiener 
filtering, Equation [T3l becomes 

X 2 = (A/V /2 7° - *(„)<!(„) ^(A/^V - *(»)<*(»)) 
+ a a U C a(n) a (n) (16) 

where C ( n ) = (a^a, n ^) and a is a tuning parameter 
which lies in the range < a < 1. Note that for a = 0, 
the result is the same as in the unfiltered case. Mini- 
mizing Equation [16] with respect to a gives the filtered 
estimator 

where we have defined M(„ jQ ) = [lift jW"\P( n ) +aA^], 
and A(„) is the truncated diagonal matrix of eigenvalues 
associated with *( n ). 

3. TESTING KL RECONSTRUCTIONS 

In this section we show results of the KL analysis of 
shear fields for a sample geometry. In Section 13.11 we 
discuss the general properties of shear KL modes for un- 
masked fields, while in Section 13.21 we discuss KL shear 
reconstruction in the presence of masking. 



3.1. KL Decomposition of a Single Field 

To demonstrate the KL decomposition of a shear field, 
we assume a square field of size 1° x 1°, divided into 
64 x 64 pixels. We assume a source galaxy density of 
20 arcmin -2 - appropriate for a ground-based survey 
such as DES - and calculate the KL basis following the 
method outlined in Section 12.11 For the computation 
of the nonlinear matter power spectrum, we assume a 
flat ACDM cosmology with Qm = 0.27 at the present 
day, with the power spectrum normalization given by 
cr 8 '= 0.81. 

Figure [1] shows a selection of nine of the 4096 shear 
eigenmodes within this framework. The KL modes 
are reminiscent of 2D Fourier modes, with higher-order 
modes probing progressively smaller length scales. This 
characteristic length scale of the eigenmodes can be seen 
quantitatively in Figure [2] Here we have computed the 
rotationally averaged power spectrum Cg for each indi- 
vidual Fourier mode, and plotted the power vertically 
as a density plot for each mode number. Because the 
KL modes are not precisely equivalent to the 2D Fourier 
modes, each contains power at a range of values in i. 
But the overall trend is clear: larger modes probe smaller 
length scales, and the modes are very close to Fourier in 
nature. 

As noted in Section [31 one useful quality of a KL de- 
composition is its diagonalization of the signal and noise 
of the problem. To explore this property, we plot in the 
upper panel of Figure [3] the eigenvalue profile of these 
KL modes. By construction, higher-order modes have 
smaller KL eigenvalues. What is more, because the noise 
in the covariance matrix is whitened (sec Scction [2.2[) . the 
expectation of the noise covariance within each mode is 
equal to 1. Subtracting this noise from each eigenvalue 
gives the expectation value of the signal-to-noise ratio: 
thus we see that the expected signal-to-noise ratio of the 
eigenmodes is above unity only for the first 17 of the 4096 
modes. 

At first glance, this may seem to imply that only the 
first 17 or so modes are useful in a reconstruction. On 
the contrary: as seen in the lower panel of Figure El these 
first 17 modes contain only a small fraction of the total 
information in the shear field (This is not an unexpected 
result: cosmic shear measurements have notoriously low 
signal-to-noise ratios!) About 900 modes are needed to 
preserve an average of 70% of the total signal, and at this 
level, each additional mode has a signal-to-noise ratio of 
below 1/10. The noisy input shear field can be exactly 
recovered by using all 4096 modes: in this case, though, 
the final few modes contribute two orders-of-magnitude 
more noise than signal. 

3.2. Testing KL Interpolation 

To test this KL interpolation technique, we use simu- 
lated shear catalog^. These catalogs contain 220 square 
degrees of simulated shear maps, computed using a ray- 
tracing grid through a cosmological N-body simulation 
of the standard A-CDM model. The shear signal is com- 
puted at the locations of background galaxies with a me- 
dian redshift of about 0.7. Galaxies are incorporated in 

3 The simulated shear catalogs were kindly made available to us 
by R. Wechslcr, M Busha, and M. Becker. 
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Fig. 1. — A sample of nine of the 4096 KL eigenmodes of a 1° X 1° patch of the sky partitioned into 64 X 64 pixels. Black is pos itive , 
red is negative, and each mode has unit norm. The modes are calculated from the theoretical shear correlation function (see Section l2.3H . 
As a conseq uenc e of the isotropy of the cosmic shear field, the covariance matrix - and thus the associated eigenmodes - are purely real 
(sec Section |3.1| I. 



the si mulation using the ADDGALS algorithm (jWechslerl 
12004 Wechsler et al. in preperation), tuned to the ex- 
pected observational characteristics of the DES mission. 

We pixelize this shear field using the same pixel size 
as above: 64 x 64 pixels per square degree. To perform 
the KL procedure on the full field with this angular res- 
olution would lead to a data vector containing over 10 6 
elements, and an associated covariance matrix contain- 
ing 10 12 entries. A full eigenvalue decomposition of such 
a matrix is computationally infcasible, so we reconstruct 
the field in 1° x 1° tiles, each 64 x 64 pixels in size. To 
reduce edge effects between these tiles, we use only the 
central 0.5° x 0.5° region of each, so that covering the 
300 square degree field requires 1200 tiles. 

In order to generate a realistic ma sk over the fiel d area , 
we follow the procedure outlined in IHikage e t al.l (|2010J l 
which generates pixel-level masks characteristic of point- 
sources, saturation spikes, and bad CCD regions. Wc 



tune the mask so that 20% of the shear pixels have no 
data. The geometry of the mask over a representative 
patch of the field can be seen in the lower panels of Fig- 
ure [U where we also show the result of the KL interpo- 
lation using 900 out of 4096 modes, with a = 0.15 (For a 
discussion of these parameter choices, see Appendix |A"| . 

The upper panels of Figure [4] give a qualitative view of 
the difficulty of cosmic shear measurements. The top left 
panel shows the noiseless shear across the field, while the 
top right panel shows the shear with shape noise for a 
DES-type survey (cr e = 0.3, n ga i = 20 arcmin -2 ). To the 
eye, the signal seems entirely washed out by the noise. 
Nevertheless, the shear signal is there, and can be fairly 
well-recovered using the first 900 KL modes (middle-left 
panel). For masked data, we must resort to the tech- 
niques of Section 12.51 to fill-in the missing data. The 
middle-right panel shows this reconstruction, with gray 
shaded regions representing the masked area. A visual 
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masses and spatial dist ribution of galaxy clusters (e.g. 
iPress fc Schechtenll974f ). Galaxy clusters are the largest 
gravitationally bound objects in the universe, and as such 
are exponentially sensitive to cosmological parameters 
(| White et al .11 19931 ). The spatial distribution of clusters 
and redshift evolution of their abundance and cluster- 
ing is sensitive to both geometrical effects of cosmol- 
ogy, as well as growth of structure. Because of this, 
cluster catalogs can be used to derive constraints on 
many interesting cosmological quantities, including the 
ma tter density Hm and power spectrum normalization 
as (|Lin et al.l 120031). the density and possible evolution 



Fig. 2. — The normalized power spectrum of each KL mode. For constant mode number, the figure represents a histogram of the power 
in that KL mode, normalized to a constant total power. KL modes represent a linear combination of Fourier modes, so that the power in 
each KL mode is spread over a range of I values. Nevertheless, the general trend is clear: larger mode numbers are associated with larger 
wave numbers, and thus smaller length scales. 

comparison of the masked and unmasked n = 900 panels 
of Figure 2] confirms qualitatively that the KL interpola- 
tion is performing as desired. This is especially apparent 
near the large cluster located at (RA,DEC)=(11.9,36.7). 
The remaining two lower panels of Figure 2] show cases 
of over-fitting and under-fitting of the shear data. If too 
few KL modes are used, the structure of the input shear 
field is lost. If too many KL modes are used, the masked 
regions are over-fit, causing the interpolated shear values 
to become unnaturally large. This observation suggests 
one rubric by which the ideal number of modes can be 
chosen; see the discussion in Appendix [A"l 

It is interesting to explore the limits of this interpo- 
lation algorithm. Figure [5] shows the KL reconstruction 
with increasing masked fractions, using a noise level typ- 
ical of space-based lensing surveys (n ga i = 100arcmin~ 2 ) 
Though the quality of the reconstruction understandably 
degrades, the lower panels show that large features can 
be recovered even with up to 50% of the pixels masked. 

For a quantitative analysis of the effectiveness of the 
KL interpolation in convergence mapping, and the poten- 
tial biases it introduces, a large-scale statistical measure 
is most appropriate. In the following sections, we test 
the utility of this KL interpolation scheme within the 
framework of shear peak statistics. 

4. SHEAR PEAK STATISTICS 

It has long been recognized that much useful 
cosmological information can be deduced from the 



of da rk energy (|Linder fc Jenkins] 12003: Vik hlinin et al 



2009 ) . primordial non -gaussianities (|Matarrese et al 



200C : iGrossi et ani2007l). and the b aryon mass fraction 
dLin et al .1 120031 : IGiodini et al.H2009l) . 

Various methods have been developed to measure the 
mass and spatial distribution of galaxy clusters, and each 
are subject to their own difficult astrophysical and ob- 
servational biases. They fall into four broad categories: 
optical or infrared richness, X-ray luminosity and sur- 
face brightness, Sunyaev-Zeldovich decrement, and weak 
lensing shear. 

While it was long thought that weak gravitational 
lensing studies would lead to robust, purely mass- 
selected cluster surveys, it has since become clear that 
shape noise and projection effects limit the usefulness 
of weak lensing in determining the 3D cluster mass 
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Fig. 3. — The eigenvalues associated with the eigenmodes discussed in Figure [T] By construction, the eigenvalue is proportional to 
the sum of signal and noise within each mode. The upper figure shows the value per mode, while the lower figure shows the normalized 
cumulative value. The stepped-pattern evident in the upper panel is due to the presence of degenerate eigenmodes which have identical 
eigenvalues (e.g. modes n = 2 a nd n = 3, related by parity as evident in Figurefl}. Because the eigenmodes arc computed from a whitened 
covariance matrix (see Section 12. 21 , the noise contribution within each mode is equal to 1. Subtracting this contribution leads to the 
plot of signal only: this shows that the signal-to-noise ratio is above unity only for the first 17 modes. Still, as seen in the lower panel, 
higher modes are required: the first 17 modes account for only 12% of the signal on average. To recover 70% of the signal in a particular 
reconstruction requires about 1000 modes. At this point, each additional mode has a signal-to-noise ratio of less than 0.1. Such a small 
signal-to-noise ratio is a well-known aspect of cosmic shear studies. 



function (lHamana et al.l)2004t iHennawi fe Spergell 120051 : 
IMandelbaum et all 120101: IVanderPlas et all 1201 lC The 

shear observed in weak lensing is non-locally related to 
the convergence, a measure of projected mass along the 
line of sight. The difficulty in deconvolving the corre- 
lated and uncorrclated projections in this quantity leads 
to difficulties in relating these projected peak heights to 
the masses of the underlying clusters in three dimen- 
sions. Recent work has shown, however, that this diffi- 
culty in relating the observed quantity to theory may be 
overcome through the use of statistics of the projected 
d ensity itself. 

iMarian et all ()2009L 120101 ) first explored the extent to 
which 2D projections of the 3D mass field trace cosmol- 
ogy. They found, rather surprisingly, that the statistics 
of the projected peaks closely trace the statistics of the 
3D peak distribution: in N - body simulations, both scale 
with the iSheth fe Tormenl (|1999f > analytic scaling rela- 
tions. The same correlated projections which bias cluster 
mass estimates contribute to a usable signal: statistics 
of projected mass alone can provide useful cosmological 
constraints, without the need for bias-prone conversions 
from peak height to cluster mass. 

A host of other work has explored diverse aspects of 
these shear peak statistics, including tests of these meth- 



ods with ensembles of N-bod y simulations (IWang et all 
l2009t iKratochvil eTaD l2010t [Dietrich fc Hartlapl 120101 )" 
the performance of various filtering functions and peak 
detec t ion statistics dPires et all l2009t iSchmidt fc Rozol 
120101 : IKratochvil et all 120111) . exploration of the spa- 
tial co rrelation of nois e with signal within convergence 
maps (|Fan et all 120101 ) , and exploration of shear-peak 
const raints on primordial non-gaussianity (jMaturi et all 
120111) . The literature has yet to converge on the ideal 
mapping procedure: convergence maps, gaussian filters, 
various matched filters, wavelet transforms, and more 
novel filters arc explored within the above references. 
There is also variation in how a "peak" is defined: sim- 
ple local maxima, "up-crossing" criteria, fractional areas 
above a certain threshold, connected-component label- 
ing, hierarchical methods, and Minkowski functionals are 
all shown to be useful. Despite diverse methodologies, all 
the above work confirms that there is useful cosmolog- 
ical information within the projected peak distribution 
of cosmic shear fields, and that this information adds to 
that obtained from 2-point statistics alone. 

4.1. Aperture Mass Peaks 

Based on this consensus, we use shear peak statistics to 
explore the possible bias induced by the KL interpolation 
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Fig. 4. — This figure illustrates the reconstruction of a small patch of masked shear from simulated shear catalog, upper panels: 
The underlying noiseless shear signal (left), the observed, noisy shear signal (middle), and the unmasked reconstruction with 900 modes 
and a = 0.15. The amplitude of the noise is calculated using an intrinsic ellipticity a e = 0.3, with an average number density of 
n gal = 20 arcmin -2 . The large peak in the upper portion of the figure is well-recovered by the KL reconstruction, lower panels: The 
KL reconstruction of the shear in the presence of 20% masking, with increasing number of modes n. The mask is rep rese nted by the 
shaded regions in panels: within these regions, the value of the shear is recovered through KL interpolation (see Section l2.5j l. We see in 
this progression the effect of the KL cutoff choice: using too few modes leads to loss of information, while using too many modes leads to 
over- fitting within the masked regions (Sec Appendix [A] for a discussion of the choice of number of modes). 



method outlin ed above. We follow the aperture mass 
methodology of lDietrich fc Hartlapl (j2010f ): The aperture 

mass magnitude at a point 9q is given by 

M ap (0o) = / d 2 6Q NFW (# = \6- 0o|)7t(0; 0o) (18) 
Jn 

where 7t(#; 9o) is the component of the shear at location 

8 tangential to the line 9 — 6q, a nd Qnfw is the NFW - 
matched filter function defined in lSchirmer et all (|2007f ): 

1 tanh(a:/x c ) 
Qnfw(x;x c ) oc - - e6 _ 150x + e _ 47+50x -^r 

(19) 

with x = i?/i>m M and x r _ a free parameter. We follow 
IDietrich fe Hartlad (|2010D and set i9 max = 5.6 arcmin 
and x c — 0.15. The integral in Equation [18] is over the 
whole sky, though the filter function Q effectively cuts 
this off at a radius $ max - In the case of our pixelized 
shear field, the integral is converted to a discrete sum 
over all pixels, with $ equal to the distance between the 
pixel centers: 

M ap (£) = ^q 

NFw($ij)7t(^; ®j) (20) 

3 

where we have defined ??y = \6i — 8j\. 



We can similarly compute the B-mode aperture 
mass, by substituting j t ~> 7x in Equations [TM2T)1 
(jCrittenden et all 12002} ). For pure gravitational weak 
lensing with an unbiased shear estimator, the B-mode 
signal is expected to be negligible, though second-order 
effects such as source clustering and intrinsic align- 
ments can cause contamination on small angu lar scales 
(jCrittenden et a,l.ll2002t ISchneider et al.ll2002b[) These ef- 
fects aside, the B-mode signal can be used as a rough 
estimate of the systematic bias of a particular analysis 
method. 

For our study, the aperture mass is calculated with 
the same resolution as the shear pixclization: 64 2 pixels 
per square degree. A pixel is defined to be a peak if its 
value is larger than that of the surrounding eight pixels: 
a simple local maximum criterion. 

4.2. The Effects of Masking 

When a shear peak statistic is computed across a field 
with masked regions, the masking leads to a bias in the 
peak height distribution (see Section[5]bclow). Moreover, 
due to the non-local form of the aperture mass statistic, a 
very large region is affected: in our case, a single masked 
pixel biases the aperture mass measurement of an area 
of size 7r(2$ max ) 2 w 400 arcmin 2 . There are two naive 
approaches one could use when measuring the aperture 
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Fig. 5. — Here wc show the same field as in Figure[4] reconstructed using n = 900 modes, with increasing levels of mask coverage. The 
density of source galaxies has been increased to 100 arcmin -2 , typical of a space-based weak lensing survey. At this noise level, smaller 
halos can be detected within the unmasked KL reconstruction (upper-left panel). Even at a 50% masking level, the large peak at (RA,DEC) 
= (11.9,36.75) is adequately recovered. 



mass in this situation: 

Unweighted:: Here we simply set the shear value within 
each masked pixel to zero, and apply Equation 1201 
The shear within the masked regions do not con- 
tribute to the peaks, so the height of the peaks will 
be underestimated. 

Weighted:: Here we implement a weighting scheme 
which re-normalizes the filter Qnfw to reflect the 
reduced contribution from masked pixels. The in- 
tegral in equation [T8l is replaced by the normalized 
sum: 

where w(6j) = if the pixel is masked, and 1 other- 
wise. This should correct for the underestimation 



of peak heights seen in the unweighted case. 

In order to facilitate comparison between this weighted 
definition of M ap and the normal definition used in the 
unmasked and unweighted cases, we normalize the lat- 
ter by QNFw(^ij); which is a constant normalization 
across the field. 

Note that in both cases, it is the shear that is masked, 
not the M ap peaks. Aperture mass is a non-local mea- 
sure, so that the value can be recovered even within the 
masked region. This means that masking will have a 
greater effect on the observed magnitude of the peaks 
than it will have on the count. In particular, on the 
small end of the peak distribution, where the peaks are 
dominated by shape noise, the masking of the shear sig- 
nal is likely to have little effect on the distribution of 
peak counts. This can be seen in Figure [71 
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Fig. 6. — Comparison of the masked and unmasked peak distributions, left panel: the peak distributions without the use of KL. The black 
line is the result with no masking, while the red and green lines show the two naive methods of correcting for the mask (see Section 14.21 . 
right panel: the masked and unmasked peak distributions after applying KL. Neither naive method of mask-correction adequately recovers 
the underlying peak distribution. It is evident, however, that the KL-based interpolation procedure recovers a mass map with a similar 
peak distribution to the unmasked KL map. It should be noted that the unmasked peak distribution (black line, left panel)is not identical 
to the unmasked peak distribution after application of KL (black line, right panel). This difference is addressed in Figure[7] 



4.3. M a p Signal-to- Noise 

It is common in shear peak studies to study signal- 
to-noise peaks rather than directl y study aperture- 
mass or convergence peaks (e.g. iWang et al.l 120091 : 
IDietrich fc Hartlad [2010l : ISchmidt fc Rozol 120101) . We 
follow this precedent here. The aperture mass (Eqn. [20)) 
is defined in terms of the tangential shear. Because 
we assume that the shear measurement is dominated by 
isotropic, uncorrclated shape noise, the noise covariancc 
of M ap can be expressed 



[Ar A/ ] lJ = (M ap (^)M ap (^)) 



(22) 



where we have used the fact that shape noise is uncorre- 
lated: [Af-y]ij = (nifij) oc 5ij. 

In the case of a KL-reconstruction of a masked shear 
field, the reconstructed shear has non- negligible correla- 
tion of noise between pixels. From Equations [T51IT71 it 
can be shown that 

= A/'y 2 *M" 1 * t VF 2 *M" 1 * t A/'y 2 (23) 

The covariance matrix Af^ is no longer diagonal, but 
the noise remains isotropic under the linear transforma- 
tion, so that A/^ t = The aperture mass noise 
covariance can thus be calculated in a similar way to the 
non-KL case: 

Wm]ij = 77 ^ ^ QNFw(^fc)QNFw(tV)[A/%]fc£. 

k i 

(24) 

This expression can be computed through standard lin- 
ear algebraic techniques. The aperture mass signal-to- 
noisc in each pixel is given by 



[S/N]i = A/ ap (0i)/vW 



M\ii 



(25) 



Unmasked Peak Distributions 
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Fig. 7. — Comparison of the distribution of M ap peaks for un- 
masked shear, before and after filtering the field with KL (dot- 
ted line and dashed line, respectively). The peak distribution in 
the absence of noise is shown for comparison (solid line). It is 
clear that the addition of shape noise leads to many spurious M ap 
peaks: noise peaks outnumber true peaks by nearly a factor of 
10 for smaller peak heights. Filtering by KL reduces these spuri- 
ous peaks by about a factor of 3, and for larger peaks leads to a 
distribution similar in scale to that of the noiseless peaks. 



5.1. 



5. DISCUSSION 
M ap Peak Distributions 



In Figures [6][9] we compare the peak distribution ob- 
tained with and without KL. We make three broad qual- 
itative observations which point to the efficacy of KL in 
interpolation of masked shear fields, and in the filtration 
of shape-noise from these fields. We stress the qualitative 
nature of these results: quantifying these observations 
in a statistically rigorous way would require shear fields 
from an ensemble of cosmology simulations, which is be- 
yond the scope of this work. These results nevertheless 
point to the efficacy of KL analysis in this context. 

KL filtration corrects for the bias due to mask- 
ing. Figure [5] compares the effect of masking on the 
resulting peak distributions with and without KL. The 
left panel shows the unmasked noisy peak distribution, 
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B-modes vs noise-only (unmasked) 




Fig. 8. — The comparison between B-mode peak distributions 
and the peak distributions for a shear field composed entirely of 
noise. As expected, the B-mode peak distributions arc largely con- 
sistent with being due to noise only. Because of this, we can use 
B-mode peaks as a rough proxy for the noise. 
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Fig. 9. — top panel: The cumulative distributions of peaks in 
signal- to- n oise, for peaks with M^/c tm > 3.25. This is the statis- 
tic used by Dietrich & Hartlap (2010) to discriminate between cos- 
mological models, bottom panel: The ratio of B-mode to E-mode 
peak distributions. Filtration by KL reduces the relative number 
of B-mode peaks by about 1/3. Because B-mode peaks are a proxy 
for contamination by shape noise (see Figurc[8]l, this indicates that 
KL-filtration results in peak distributions less affected by statisti- 
cal errors. KL also reduces the total number of both E and B peaks 
by about 2/3; this effect can also be seen in Figure 

and the masked peak distributions resulting from the 
weighted and unweighted approaches described in Sec- 
tion 14.21 Neither method of accounting for the masking 
accurately recovers the unmasked distribution of peak 
heights. The unweighted approach (green line) leads to 
an underestimation of peak heights. This is to be ex- 
pected, because it does not correct for the missing infor- 
mation in the masked pixels. The weighted approach, on 
the other hand, over-estimates the counts of the peaks. 
We suspect this is due to an analog of Eddington bias: 
the lower signal-to-noise ratio of the weighted peak statis- 
tic leads to a larger scatter in peak heights. Because 
of the steep slope of the peak distribution, this scatter 
preferentially increases the counts of larger peaks. This 
suspicion is confirmed by artificially increasing the noise 
in the unmasked peak function. Increasing er e from 0.30 
to 0.35 in the unmasked case results in a nearly identical 
peak function to the weighted masked case. 
The right panel of Figure [6] shows that when KL is ap- 



plied to the shear field, the distribution of the masked 
and unmasked peaks is very close, both for E-mode and 
B-mode peaks. This indicates the success of the KL- 
based interpolation outlined in Section 12.51 Even with 
20% of the pixels masked, the procedure can recover 
a nearly identical peak distribution as from unmasked 
shear. 

KL filtration reduces the number of noise peaks. 

Comparison of the unmasked lines in the left and right 
panels of Figure[6]shows that application of KL to a shear 
field results in fewer peaks at all heights. This is to be 
expected: when a reconstruction is performed with fewer 
than the total number of KL modes, information of high 
spatial frequency is lost. In this way, KL acts as a sort 
of low-pass filter tuned to the particular signal-to-noise 
characteristics of the data. Figure [7] over-plots the KL 
and non-KL peak distributions with the noiseless peak 
distribution. From this figure we see that the inclusion of 
shape noise results in nearly an order-of-magnitude more 
peaks than the noiseless case. The effect of noise on peak 
counts lessens slightl y for higher-M an p e aks: this sup- 
ports the decision of iDietrich fc Hartlapl (|2010f) to limit 
their distributions to peaks with a signal-to-noise ratio 
greater than 3.25: the vast majority of peaks are lower 
magnitude, and are overwhelmed by the effect of shape 
noise. 

Omission of higher-order KL modes of shear field re- 
duces the number of these spurious peaks by a factor of 
3 or more. For low- magnitude peaks, (M ap < 0.02), KL 
still produces peak counts which are dominated by noise. 
For higher-magnitude peaks, the number of observed KL 
peaks more closely approaches the number of peaks in 
the noiseless case. 

KL filtration reduces the presence of B-modes. 
To first order, weak lensing shear is expected to consist 
primarily of curl-free, E-mode signal. Because of this, 
the presence of B-modes can indicate a systematic effect. 
It is not obvious that filtration by KL will maintain this 
property: as noted in Scction [2T4l KL modes individually 
are agnostic to E-mode and B-modc information. E&B 
information is only recovered within a complex-valued 
linear combination of the set KL modes. 

Figure |S] shows a comparison between the unmasked 
B-mode peak functions from Figure [S] and the associ- 
ated E-mode peak functions due to shape-noise only. For 
both the non-KL version and the KL version, the B-mode 
peak distributions closely follow the distributions of noise 
peaks. This supports the use of B-mode peaks as a proxy 
for the peaks due to shape noise, even when truncating 
higher-order KL modes. 

The near-equivalence of B-modes and noise-only peaks 
shown in Figure [8] suggests a way of recovering the true 
peak function, by subtracting the B-mode count from the 
E-mode count as a proxy for the shape noise. This ap- 
proach has one fatal flaw: because it involves computing 
the small difference between two large quantities, the re- 
sult has extremely large uncertainties. It should be noted 
that this noise contamination of small peaks is not an im- 
pediment to using this method for cosmological analyses: 
the primary information in shear peak statistics is due 
to the high signal-to-noise peaks. 

In the top panel of Figure |H1 we show the cumulative 
distribution of peaks in signal-to-noise, for peaks with 
S/N > 3.25: the quantity used as a cosmological dis- 
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criminant in iDietrich fc Hartlad pOlOh . The difference 
in the total number of E-mode peaks in the KL and 
non-KL approaches echoes the result seen in Figure [7j 
truncation of higher-order KL modes acts as a low-pass 
filter, reducing the total number of peaks by a factor of 
~ 3. More interesting is the result shown in the lower 
panel of Figure [SJ where the ratio of B-mode peak counts 
to E-mode peak counts is shown. Before application of 
KL, the B-mode contamination is above 30%. Filtra- 
tion by KL reduces this contamination by a factor of 
~ 3, to about 10%. This indicates that the truncation of 
higher-order KL modes leads to a preferential reduction 
of the B-mode signal, which traces the noise. This is a 
promising observation: the counts of high signal-to-noise 
peaks, whic h offer the most sensitivit y to cosmological 
parameters ([Dietrich fc Hartrapll2010f ). are significantly 
less contaminated by noise after filtering and reconstruc- 
tion with KL. This is a strong indication that the use of 
KL could improve the cosmological constraints derived 
from studies of shear peak statistics. 

Note that in Figure [9] we omit the masked results 
for clarity. The masked cumulative signal-to-noise peak 
functions have B/E ratios comparable to the unmasked 
versions, so the conclusions here hold in both the masked 
and unmasked cases. 



5.2. Remaining Questions 

The above discussion suggests that KL analysis of 
masked shear fields holds promise in constraining cos- 
mological parameters of shear peaks in both masked and 
unmasked fields. KL greatly reduces the number of spuri- 
ous noise peaks at all signal-to-noise levels. It minimizes 
the bias between masked and unmasked constructions, 
and leads to a factor of 3 suppression of the B-mode sig- 
nal, which is a proxy for the spurious signal introduced 



through shape noise. 

The question remains, however, how much cosmolog- 
ical information is contained in the KL peak functions. 
The reduction in level of noise peaks is promising, but the 
omission of higher-order modes in the KL reconstruction 
leads to a smoothing of the shear field on scales smaller 
than the cutoff mode. This smoothing could lead to the 
loss of cosmologically useful information. In this way, 
the choice of KL mode cutoff can be thought of as a bal- 
ance between statistical and systematic error. The effect 
of these competing properties on cosmological param- 
eter determination is difficult to estimate. Quantifying 
this effect will require analysis within a suite of synthetic 
shear maps, similar to the approach t aken in previous 
studi es (e.g. IDietrich fc Hartlaa 120101 : iKratochvil et al.l 
I2010I ). and will be the subject of future work. 

Another possible application of KL in weak lensing is 
to use KL to directly constrain 2-point information in 
the measured shear data. In contrast to the method out- 
lined in the current work, KL basis functions can be com- 
puted for the unmasked region only. The projection of 
observed data onto this basis can be used to directly com- 
pute cosmological parameters via the 2-point function, 
without ever explicitly calculating the power spectrum. 
T his is similar to the appr oach taken for galaxy counts 
in lVogelev fc Szalavl (|1996( ). We explore this approach in 
a followup work. 
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APPENDIX 
CHOICE OF KL PARAMETERS 

The KL analysis outlined in Section[2]has only two free parameters: the number of modes n and the Wiener filtering 
level a. Each of these parameters involves a trade-off: using more modes increases the amount of information used 
in the reconstruction, but at the expense of a decreased signal-to-noise ratio. Decreasing the value of a to reduces 
the smoothing effect of the prior, but can lead to a nearly singular convolution matrix iVf^ Q,), which results in 
unrealistically large shear values in the poorly-constrained areas areas of the map (i.e. masked regions). 

To inform our choice of the number of modes n, we recall the trend of spatial scale with mode number seen in 
Figure [2j Our purpose in using KL is to allow interpolation in masked regions. To this end, the angular scale of the 
mask should inform the choice of angular scale of the largest mode used. An eigcnmode which probes scales much 
smaller than the size of a masked region will not contribute meaningful information to the reconstruction within that 
masked region. Considering the pixels within our mask, we find that 99.5% of masked pixels are within 2 pixels of a 
shear measurement. This corresponds to an angular scale of £ = 6140. Consulting Figure [21 we see that modes larger 
than about n = 900 out of 4096 will probe length scales significantly smaller than the mask scale. Thus, we choose 
n = 900 as an appropriate cutoff for our reconstructions. 

To inform our choice of the Wiener filtering level a, we examine the agreement between histograms of A/ ap peaks 
for a noise-only DES field with and without masking (see Section [4]). We find that for large (small) values of a, the 
number of high-M ap peaks is underestimated (overestimated) in the masked case as compared to the unmasked case. 
Empirically, we find that the two agree at a = 0.15; we choose this value for our analysis. Note that this tuning is done 
on noise-only reconstructions, which can be generated for observed data by assuming that shape noise dominates: 

Wi\a = ^ ij . (Al) 

The a-tuning can thus be performed on artificial noise realizations which match the observed survey characteristics. 

We make no claim that (n, a) = (900, 0.15) is the optimal choice of free parameters for KL: determining this would 
involve a more in-depth analysis. They are simply well-motivated choices which we use to make a case for further 
study. 



